Load network and scent data

From Tables S1 and S2 in: Kantsa, A., R.A. Raguso, A.G. Dyer, J.M. Olesen, T. Tscheulin, and T. Petanidou. 2018. Disentangling the role of floral sensory stimuli in pollination networks. Nature Communications 9: 1041. doi:[10.1038/s41467-018-03448-w](https://doi.org/10.1038/s41467-018-03448-w)

knitr::opts_chunk$set(echo = TRUE, comment="  ", cache=TRUE, fig.width=14, fig.asp=.678)
library(openxlsx)
library(data.table)
library(RColorBrewer)
pal <- brewer.pal(5,"Set2")
library(heatmaply)
## Loading required package: plotly
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
## Loading required package: viridis
## Loading required package: viridisLite
## 
## ======================
## Welcome to heatmaply version 0.14.1
## 
## Type citation('heatmaply') for how to cite the package.
## Type ?heatmaply for the main documentation.
## 
## The github page is: https://github.com/talgalili/heatmaply/
## Please submit your suggestions and bug-reports at: https://github.com/talgalili/heatmaply/issues
## Or contact: <tal.galili@gmail.com>
## ======================
library(bipartite)
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-1
## Loading required package: sna
## Loading required package: statnet.common
## 
## Attaching package: 'statnet.common'
## The following object is masked from 'package:base':
## 
##     order
## Loading required package: network
## network: Classes for Relational Data
## Version 1.13.0.1 created on 2015-08-31.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
##                     Mark S. Handcock, University of California -- Los Angeles
##                     David R. Hunter, Penn State University
##                     Martina Morris, University of Washington
##                     Skye Bender-deMoll, University of Washington
##  For citation information, type citation("network").
##  Type help("network-package") to get started.
## sna: Tools for Social Network Analysis
## Version 2.4 created on 2016-07-23.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
##  For citation information, type citation("sna").
##  Type help(package="sna") to get started.
##  This is bipartite 2.08
##  For latest changes see versionlog in  ?"bipartite-package".
##  For citation see: citation("bipartite").
##  Have a nice time plotting and analysing two-mode networks.
## 
## Attaching package: 'bipartite'
## The following object is masked from 'package:vegan':
## 
##     nullmodel
data(mosquin1967)
library(vegan)
library(phytools)
## Loading required package: ape
## 
## Attaching package: 'ape'
## The following object is masked from 'package:sna':
## 
##     consensus
## Loading required package: maps
library(taxize)

#load network
net <- read.delim("Kantsa_etal_Supp_Data1.csv", skip = 2)
rownames(net) <- net[,1]
net[,1] <- NULL
net <- as.matrix(t(net))
snet <- sortweb(net) #sort by row & column totals

#load floral scents
scent.file <- "Kantsa_2018_Supp_Data2.xlsx"
plants <- getSheetNames(scent.file)
scentlist <- lapply(2:length(plants), read.xlsx, xlsxFile=scent.file)
plants <- plants[-1]
scent <- rbindlist(scentlist, fill=T, idcol="plant")[,-4]
colnames(scent)[2:3] <- c("cmpd","mean")
scent$plant <- factor(plants[scent$plant])
scent$cmpd <- factor(tolower(scent$cmpd))
cmpds <- levels(scent$cmpd)

scent.net <- dcast(scent, plant ~ cmpd, value.var="mean", fun.aggregate=mean, fill=0) #long to wide
scent.net <- as.data.frame(scent.net)
rownames(scent.net) <- scent.net[,1]
scent.net[,c(1,380)] <- NULL #weird NA column has one entry
scent.net <- as.matrix(scent.net)

Plant phylogeny

#plants.class <- classification(plants, db="gbif")
#save(plants.class, file="plantsclass.Rdata")
load("plantsclass.Rdata")
plants.tree <- class2tree(plants.class)
plot(plants.tree)

Plant-pollinator network

Heatmap

heatmaply(net, scale="column")

Matrix plot

#interaction matrix plot
visweb(net, type="nested", circles=T, boxes=F, circle.max=1.2)

Web plot

#interaction web plot
plotweb(net, text.rot=90, method="cca", col.high=pal[2], col.low=pal[1], col.interaction=pal[5])

Network topography

networklevel(net)
                         connectance                     web asymmetry 
                          0.06312657                        0.63106796 
                   links per species            number of compartments 
                          1.95631068                        2.00000000 
               compartment diversity               cluster coefficient 
                          1.07901412                        0.02631579 
                          nestedness               weighted nestedness 
                          4.48059576                        0.58672881 
                       weighted NODF    interaction strength asymmetry 
                         10.80348739                        0.20372363 
            specialisation asymmetry                   linkage density 
                         -0.04424002                        3.57255189 
                weighted connectance                      Fisher alpha 
                          0.01734248                       89.61252715 
                   Shannon diversity              interaction evenness 
                          3.27629195                        0.37393976 
        Alatalo interaction evenness                                H2 
                          0.25285628                        0.57919894 
                number.of.species.HL              number.of.species.LL 
                        168.00000000                       38.00000000 
   mean.number.of.shared.partners.HL mean.number.of.shared.partners.LL 
                          0.34944397                        1.18776671 
              cluster.coefficient.HL            cluster.coefficient.LL 
                          0.34062967                        0.21733898 
     weighted.cluster.coefficient.HL   weighted.cluster.coefficient.LL 
                          0.72421977                        0.31277076 
                    niche.overlap.HL                  niche.overlap.LL 
                          0.12366249                        0.14889538 
                     togetherness.HL                   togetherness.LL 
                          0.03019814                        0.02302573 
                          C.score.HL                        C.score.LL 
                          0.75923950                        0.74100046 
                          V.ratio.HL                        V.ratio.LL 
                          3.23242952                       17.02718828 
                      discrepancy.HL                    discrepancy.LL 
                        248.00000000                      250.00000000 
                 extinction.slope.HL               extinction.slope.LL 
                          5.03260774                        1.60178636 
                       robustness.HL                     robustness.LL 
                          0.81123903                        0.61539429 
       functional.complementarity.HL     functional.complementarity.LL 
                       8044.31456516                     7651.23911634 
                partner.diversity.HL              partner.diversity.LL 
                          0.94103312                        1.28483306 
                       generality.HL                  vulnerability.LL 
                          2.67434342                        4.47076037

Scent network

Heatmap

heatmaply(scent.net, scale="column")

Web plot

#interaction web plot
plotweb(scent.net, text.rot=90, method="cca", col.high=pal[2], col.low=pal[1], col.interaction=pal[5])

NMDS

scent.mds <- metaMDS(sqrt(scent.net), autotransform = F)
   Run 0 stress 0.2177614 
   Run 1 stress 0.2155518 
   ... New best solution
   ... Procrustes: rmse 0.08183959  max resid 0.3861441 
   Run 2 stress 0.2277 
   Run 3 stress 0.2198936 
   Run 4 stress 0.2217743 
   Run 5 stress 0.2139972 
   ... New best solution
   ... Procrustes: rmse 0.1078209  max resid 0.3550252 
   Run 6 stress 0.2201085 
   Run 7 stress 0.2306285 
   Run 8 stress 0.2158631 
   Run 9 stress 0.2417021 
   Run 10 stress 0.2264037 
   Run 11 stress 0.2330208 
   Run 12 stress 0.2175837 
   Run 13 stress 0.2506198 
   Run 14 stress 0.2138077 
   ... New best solution
   ... Procrustes: rmse 0.07751348  max resid 0.2765495 
   Run 15 stress 0.2130163 
   ... New best solution
   ... Procrustes: rmse 0.08566332  max resid 0.3435179 
   Run 16 stress 0.2426195 
   Run 17 stress 0.2136397 
   Run 18 stress 0.2156501 
   Run 19 stress 0.2146277 
   Run 20 stress 0.2169694 
   *** No convergence -- monoMDS stopping criteria:
       20: stress ratio > sratmax
scent.op <- ordiplot(scent.mds, type="n")
points(scent.op, what="species", pch=3, col=pal[3])
text(scent.op, what="sites", cex=0.8)

Phylochemospace

scent.mds.taxa <- scent.mds$points
rownames(scent.mds.taxa) <- plants.tree$phylo$tip.label

plants.tree$phylo$edge.length <- plants.tree$phylo$edge.length+0.00001

phylomorphospace(plants.tree$phylo, scent.mds.taxa, control=list(col.node=0), label="none", lwd=1)
text(scent.mds.taxa[,1], scent.mds.taxa[,2], rownames(scent.mds.taxa), cex=0.8, offset=0.5, xpd=T)

Phylo Bar chart

par(mfrow=c(1,2))
par(mar=c(0.8,0,0.8,0))
plot(plants.tree, cex=1)
par(mar=c(0,0,0,0))
barplot(t(scent.net), col=sample(rainbow(ncol(scent.net))), names.arg=rep("",nrow(scent.net)), horiz=TRUE) 

par(mfrow=c(1,2))
par(mar=c(0.8,0,0.8,0))
plot(plants.tree, cex=1)
par(mar=c(0,0,0,0))
barplot(t(decostand(scent.net, method="tot")), col=sample(rainbow(ncol(scent.net))), names.arg=rep("",nrow(scent.net)), horiz=TRUE)